Import and visualize 3d

Now that I have plot points and treatment boundaries all defined with consistent naming, I can visualize a 3d model of our topography and determine a strategy for extracting topograpical variables and simulating representative topogrphy for fire simulation.

Code

library(dplyr)
library(sf)
library(stars)
library(rayshader)

rgl::setupKnitr(autoprint = FALSE)

# I want to overlay my plots on a 3d surface representing the terrain for each
# site
aerial_files <- list.files("../data/gis/EE_naip_jdsf", full.names = TRUE, pattern = ".tif$")
dem_files <- list.files("../data/gis/usgs_1m_dem", full.names = TRUE, pattern = ".tif$")

aerial <- st_mosaic(aerial_files) |> read_stars()
dem <- st_mosaic(dem_files) |> read_stars()

vect_db <- "../data/gis/named_numbered.gpkg"
blocks_buffer <- st_read(vect_db, "blocks_buffer", as_tibble = TRUE, quiet = TRUE)
blocks <- st_read(vect_db, "blocks", as_tibble = TRUE, quiet = TRUE)
plot_centers <- st_read(vect_db, "plot_centers", as_tibble = TRUE, quiet = TRUE)
plot_discs <- st_buffer(plot_centers, 11.28)

# try it out with one block
oneblock <- filter(blocks, site == "camp8", treatment == "ls", burn == "b")
oneplots <- filter(plot_discs, site == "camp8", treatment == "ls", burn == "b")
onecenters <- filter(plot_centers, site == "camp8", treatment == "ls", burn == "b")
oneraster <- dem[oneblock]
oneimage <- aerial[oneblock]

dem_matrix <- st_as_stars(oneraster)[[1]]

# Get actual data for 3 bands scaled to [0,1], and transposed like rayrender
# likes it. colors stretched to data range to brighten
aerial_array <- st_as_stars(oneimage) |>
  slice(band, 1:3) |>
  pull() |>
  (`/`)(255) |>
  scales::rescale(to = c(0, 1)) |>
  aperm(c(2, 1, 3))

data_plots <- oneplots |> filter(!is.na(plotid))
plots_overlay <- rayshader::generate_line_overlay(
  st_cast(data_plots, "MULTILINESTRING"),
  extent = st_bbox(oneraster),
  heightmap = dem_matrix,
  color = "red"
)
rw_texture <- create_texture("#bca5a1", "#bf754e", "#6a3b1e", "#bf754e", "#daaa75")

aerial_array |> 
  add_overlay(plots_overlay) |>
  plot_3d(dem_matrix)

rgl::rglwidget()

Code

# dem_matrix |>
#   sphere_shade(texture = rw_texture) |>
#   add_overlay(plots_overlay) |>
#   plot_3d(dem_matrix)

bb <- st_bbox(oneraster)
deltay <- bb[4] - bb[2]
y_offset <- (1 - 200 / deltay) / 2 
dem_matrix |>
  sphere_shade(texture = rw_texture) |>
  add_shadow(texture_shade(dem_matrix)) |>
  add_shadow(lamb_shade(dem_matrix), max_darken = 0.2) |>
  add_overlay(plots_overlay) |>
  plot_3d(dem_matrix)
render_compass(position = "S", altitude = min(dem_matrix, na.rm = TRUE))
render_scalebar(
  limits = c(0, 200),
  scale_length = c(y_offset, 1 - y_offset),
  y = min(dem_matrix, na.rm = TRUE),
  label_unit = "m"
)
render_camera(theta = 300, phi = 25, fov = 60, zoom = 0.8)

data_centers <- onecenters |> filter(!is.na(plotid))
for (i in seq_len(nrow(data_centers))) {
  render_label(
    dem_matrix,
    text = as.character(data_centers[i, "plotnum", drop = TRUE]),
    lat = st_coordinates(data_centers)[i, 2],
    long = st_coordinates(data_centers)[i, 1],
    altitude = 10,
    extent = st_bbox(oneraster),
    linewidth = 1
  )
}

rgl::rglwidget()